Optimal control of number squeezing in trapped Bose-Einstein condensates 



Julian Grond/ Gregory von Winckel,^ Jorg Schmiedmayer,"^ and Ulrich Hohenester^ 

^ Institut fiir Physik, Karl-Franzens-Universitdt Graz, Universitdtsplatz 5, 8010 Graz, Austria 
^Institut fiir Mathematik und Wissenschaftliches Rechnen, 
Karl-Franzens-Universitdt Graz, Heinrichstrafie 36, 8010 Graz, Austria 
^ Atominstitut der Osterreichischen Universitdten, TU-Wien, Stadionallee 2, 1020 Wien, Austria 

(Dated: December 2, 2009) 

We theoretically analyze atom interferometry based on trapped ultracold atoms, and employ 
optimal control theory in order to optimize number squeezing and condensate trapping. In our 
simulations, we consider a setup where the confinement potential is transformed from a single to 
a double well, which allows to split the condensate. To avoid in the ensuing phase-accumulation 
stage of the interferometer dephasing due to the nonlinear atom-atom interactions, the atom number 
fluctuations between the two wells should be sufficiently low. We show that low number fluctuations 
(high number squeezing) can be obtained by optimized splitting protocols. Two types of solutions are 
found: in the Josephson regime we flnd an oscillatory tunnel control and a parametric ampliflcation 
of number squeezing, while in the Fock regime squeezing is obtained solely due to the nonlinear 
coupling, which is transformed to number squeezing by peaked tunnel pulses. We study splitting 
and squeezing within the frameworks of a generic two-mode model, which allows us to study the 
basic physical mechanisms, and the multi-configurational time dependent Hartree for bosons method, 
which allows for a microscopic modeling of the splitting dynamics in realistic experiments. Both 
models give similar results, thus highlighting the general nature of these two solution schemes. We 
flnally analyze our results in the context of atom interferometry. 



PACS numbers: 03.75.-b,02.60.Pn,37.25.+k,81.16.Ta 

I. INTRODUCTION 

Squeezed states offer the possibility of measurements 
below shot noise [1 , 2 . In matter wave interferometry [3^ 
with Bose Einstein condensates (BECs), number squeez- 
ing in the split BEG helps to slow down the phase dif- 
fusion arising from atom-atom interactions [J, and thus 
the potential measurement time and with it the sensitiv- 
ity of the atom interferometer is increased. 

Splitting a trapped BEG is achieved by transforming 
the confinement potential from a single to a double well 
O H]. Glose to the sphtting point, where the conden- 
sate breaks apart into two spatially separated parts, the 
competition between tunneling and nonlinear interaction 
leads to a decrease of number fluctuations and squeezing. 

In this paper we will discuss our concepts to optimize 
squeezing and interference in trapped BEGs in the con- 
text of Atom chips E [9] which provide a versatile tool 
for precise manipulation of trapped BEGs and allow ro- 
bust implementation of the splitting procedure and inter- 
ference experiments [6l [10] . A key ingredient in our con- 
sideration is the beam splitter for the trapped or guided 
atoms pT| [T2| [T3| [T4] . As an example we will consider 
here the RF beam splitter (TH [TSl [iB] as employed in the 
atom chip experiments. However, our considerations can 
be easily apphed to any other double well splitting and 
trapped atom interference setup [5| (THl [HI [TH] . 

Squeezing in the context of BEG splitting has been 
frequently discussed by assuming an approximately ex- 
ponential decrease of the tunnel couphng close to the 
sphtting point [T9( |20l [211 [22]. In experiment, number 
squeezing has been achieved in an optical trap in 



optical lattices [HI [HI [IHl [23 [SB], and more recently 
through hnear splitting [18l [29]. In Ref. [30] we have 
recently shown that exponential splitting is by far not 
optimal, and there exist alternative squeezing protocols 
allowing for much more squeezing. These solutions are 
obtained within the framework of optimal control theory 
(OGT) [31] [32]. In this paper we give details on our work 
about optimal number squeezing. 

In our simulations we go beyond the simple two-mode 
model and consider additionally the spatial dynamics, 
in order to have a reahstic description of the splitting 
dynamics and to evaluate the importance of condensate 
excitations. Our theoretical approach is based on the 
multi-configurational time- dependent Hartree for bosons 
(MGTDHB) method [32], which allows to describe both 
atom number and spatial dynamics in a self-consistent 
fashion. This method allows us to use as a control pa- 
rameter the splitting distance of the trap used in exper- 
iments. The generic two-mode model, whose control is 
through the tunnel coupling which is determined indi- 
rectly through the conflnement potential, is used to ob- 
tain a simpler interpretation of the underlying control 
strategies. 

We flnd that there exist even simpler splitting proto- 
cols compared to our earlier work [30 . These solutions 
can be understood similarly to the control schemes of 
Refs. [341 [35l |36l [37] which rely on a constant, optimal 
value of the tunnel coupling. We compare the OGT re- 
sults to a two-parameter optimization, where the control 
consists of an initial exponential splitting followed by a 
constant value of the tunnel couphng. While within the 
generic two-mode model the results almost coincide, this 
simple strategy is not generally applicable within a real- 
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istic model for two main reasons: first, for short splitting 
times (on a time scale determined by the trapping poten- 
tial) the condensates oscillate, and the split condensate 
ends up in an excited state; second, after number squeez- 
ing the condensates have to be further separated to turn 
off any tunnel coupling, which for realistic traps can only 
be achieved with refined protocols. To avoid these prob- 
lems, we are seeking for control protocols where at the 
final time the condensates are at rest and fully decoupled. 

We have organized our paper as follows. In Sec. II we 
give a brief description of our model system and its dy- 
namic equations. Squeezing and squeezing optimization 
within a generic two-mode model is described in Sec. III. 
In Sec. IV we introduce our approach for the realistic 
simulation of condensate splitting, and discuss splitting 
and its optimizations within the MCTDHB framework. 
We show that the general trends of the two-mode model 
and MCTDHB are in agreement. Finally, in Sec. V we 
analyze our optimized splitting procedures in the context 
of atom interferometry. 



II. CONDENSATE SPLITTING 
A. Model 

In our theoretical approach, we consider an elongated 
condensate whose wave function is modified only along a 
single spatial direction x [M]. The atom dynamics can 
be described by the Hamiltonian [38, 39^ 



H = 



dx . 



(1) 



The first term on the right-hand side (rhs), which com- 
prises the single particle or bare Hamiltonian h{x) = 
— ^V^ Vx{x), accounts for the kinetic energy and the 
magnetic confinement potential Vx{x). The control pa- 
rameter X(t) describes the variation of the confining po- 
tential when changing the external parameters, such as 
currents or rf-fields [8, 14 . Through X{t) it is possible 
to manipulate the trapped Bose-Einstein condensate, e.g. 
to split it by varying the potential from a single to a dou- 
ble well, see Fig. [T] The second term on the rhs accounts 
for the atom- atom interactions, which is scaled with the 
ID-interaction parameter Uq. The field operators ^(x) 
and "^^x) obey the usual equal-time commutation rela- 
tions. 

We consider ^"^Rb atoms, and use units where h = 1, 
the mass of the atoms is set to one and length is measured 
in micrometers. The natural unit of time (energy) is then 
1.37 ms (5.58 nano Kelvin), and the results for the generic 
model can be arbitrarily scaled. 



B. Confinement potential 

For the realistic modehng of the sphtting process, we 
consider the magnetic confinement potential on the chip 



modified by rf-dressing as discussed in Lesanovsky et 
al. [14^. The magnetic confinement is based on an elon- 
gated Z-wire trap, the rf field controls the splitting pro- 
cess. Throughout we use the same parameters as in 
Ref [14 , and parametrize the rf-field Bj-f = (0.5 -|- 0.3 x 
A) G through the control parameter A that determines 
the splitting distance. Note that this potential is approx- 
imated well by the function ax'^ bx"^, with A-dependent 
parameters a and b. Within this approach it is possible 
to directly determine the experimentally accessible con- 
trol field, which is often related to the tunnel coupling in 
a highly non-trivial fashion [40] . 



C. Condensate dynamics 

Since the number of atoms in the condensate is usu- 
ally large, the many-body Hamiltonian of Eq. ([T]) is ap- 
proximated by expanding the field operators in the basis 
of orbitals. A famous example is the mean-field Gross- 
Pitaevskii (GP) equation ^SSj, which is obtained by using 
a single orbital only. While this model describes well the 
dynamics of the condensate density, it does not account 
properly for quantum fluctuations. A widely used ap- 
proximation for condensates in a double- well is provided 
by the generic two- mode model [19', IT, which is ob- 
tained by restricting the field operator to a left and a 
right orbital (/)l{x) and (/)r{x), respectively. Hence, 



^(x) = dL(t)L{x) + dR(f)R{x) . 



(2) 



Within this approximation the dynamics is only related 
to the distribution of atoms between the left and right 
well, determined mainly by the tunnel coupling which is 
given as Q = — J dx(i)\(x)h(x)4>R(x) + h.c. One usually 
assumes that Q can be modified by slightly varying the 
confinement potential of the condensate, and thus Q al- 
lows to control the atom number dynamics. In contrast, 
the nonlinear intra- well energy k = ^ J dx\(l)L,Rix)\'^ is 
assumed to be not affected by such variations. This is a 
reasonable assumption since Q accounts for a tunneling 
process, whose efficiency depends extremely sensitive on 
the details of the potential, whereas k is expected to be 
rather insensitive to slight modifications of Va. 

For a more complete modeling of the whole splitting 
process, including the spatial dynamics, the ansatz for 
the field operator is 



if(x) = dLit)(f)L{x, t) + dR{t)(f)R{x, t) , 



(3) 



where the orbitals depend on time now. The number 
distribution and orbitals are calculated self-consistently 
using the multi configurational time dependent Hartree 
for Bosons method (MCTDHB) [331 • The control within 
this approach is now directly given by the trapping po- 
tential. 

In the following, we discuss optimal condensate spht- 
ting within these two models. The simpler generic model 
will provide us with a more transparent picture of the 
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FIG. 1: (Color online) Schematics of the splitting process. 
When splitting a condensate, it breaks up into two parts (red 
solid lines) with equal mean atom number, and number fluc- 
tuations An. The distribution of atoms between the wells is 
controlled by the tunnel coupling 17 and the nonlinear intra- 
well energy k. 

underlying physics, whereas the MCTDHB approach will 
allow us for a reahstic modeling of the splitting process. 



III. OPTIMIZATION OF A GENERIC 
TWO-MODE MODEL 

A. Model 

We start by discussing optimal squeezing within a 
generic two-mode model [H]. The Hamiltonian is 
obtained by plugging Eq. Q into Eq. ([T]), and reads 

^ 17 (t) /^f^ ^ ( '-i ^ .^t^t^ ^ \ 

H = —(alaii + aLa'ii) + K(a'^alaLaL + a'^a'^aRaii) . 

(4) 

In accordance to Ref. [41 , we have neglected all inter- 
well couplings, which are usually much smaller due to the 
small orbital overlap. The two- mode model is then fully 
characterized by the tunnel coupling 0(t). We choose 
units such that 2kN = 1. 

Eq. ^ can be brought into a more intriguing form 
by introducing pseudospin operators Jx, Jy, and Jz [H]. 
Jx = ^{a'^L^R + a^RttL) accounts for tunneling between 
the left and the right well, and Jz = \(a}]^aL — cl^rClr) 
measures the atom number difference between the wells. 
We label our basis states as = \N/2 + k)L\N/2-k)R, 
where is the total number of atoms and the subscripts 
denote particles in the left and right well, respectively. In 
this basis, the operators act on a state vector C whose 
entries are the probability amplitudes for a certain distri- 
bution of atoms between the wells. The two-mode Hamil- 
tonian 

H = -n{t)Jx + 2nJl (5) 

is then completely analogous to the Josephson Hamilto- 
nian for superconductors, where VL{t) is associated with 
the (time-dependent) Josephson energy and the nonhn- 
earity k, with the charging energy. 

For the unsplit trap we assume that all atoms reside in 
the bonding orbital + 4^r- Thus, the initial C vector 



FIG. 2: (Color online) Bloch sphere representation of C. 
The z-axis corresponds to the number difference (An and A0 
are number and phase fluctuations, respectively). The action 
of the pseudo-spin operators Jx and Jz on C corresponds 
to rotations about the x- and 2;-axis, respectively. In (a) a 
binomial state is shown and in (b) a number squeezed one. 

corresponds to a binomial distribution for the relative 
atom number between left and right orbital [19], ^ 




The standard deviation of the relative atom number be- 
tween left and right well associated with this state is 

Ano := {Jz) = When the trap is split, the non- 

linearity K becomes comparable to or larger than and 
number fluctuations in the ground state are reduced ow- 
ing to the cost of (nonlinear) energy. Finally, for the split 
trap we have k,:^ Q and the ground state is the number 
eigenstate | /c = 0) , which can be regarded as a completely 
squeezed state. 

A convenient representation of the atom number dis- 
tribution vector C is on the Bloch sphere [43j , Fig. [2] A 
state on the north pole would correspond to all atoms re- 
siding in the left well, and a state on the south pole to all 
atoms in the right well. All atoms in the bonding orbital 
corresponds to a state locahzed around x = 1. This bino- 
mial state has equal uncertainty in number difference {z- 
axis) and in the conjugate phase observable (around the 
equator of the sphere). In contrast, the squeezed state 
shown in panel (b) has reduced number and enhanced 
phase fluctuations. 

B. Quasi-adiabatic splitting 

Our goal is to split the condensates such that max- 
imal number squeezing is achieved. The usual strat- 



^ In this investigation, we take Q{0) for the initial state from the 
self-consistent solutions of MCHB(2) calculations 42 . This will 
allow us for a more direct comparison with the MCTDHB re- 
sults presented in Sec. |IV| Our results, however, do not depend 
significantly on the exact initial state. Due to a slight quantum 
depletion caused by interactions we obtain a correction to the 
initial squeezing. 
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FIG. 3: (Color online) Number squeezing for quasi-adiabatic 
exponential (gray, many lines), two-parameter (blue) and 
Fock optimization (symbols), for different atom numbers. For 
the optimized solutions we only report the final value of the 
number fluctuations. Both optimization strategies allow for 
number squeezing at least one order of magnitude faster than 
exponential splitting. The green lines report the limitations 
of the two-parameter optimization for = 10000 (limita- 
tions are very similar for smaller N). Solutions for realistic 
splitting within MCTDHB are only obtained for sufficiently 
long exponential decay constants tc, reported in the panel. 
Smaller values of tc result in strongly oscillating condensates 
and reduced number squeezing. We note that in all our fig- 
ures time is measured in units of 1.37 ms and position in fim 
(Rubidium units). For the generic model we use units such 
that 2kN = 1. 



egy to obtain squeezing is splitting slowly, such that the 
condensate follows adiabatically the ground state of the 
trap. The gray hues in Figure [3] show results for a quasi- 
adiabatic exponential switching-off of the tunneling rate 
according to Q{t) = O(0)e~*/*^ [19 . Depending on tc, 
at some point adiabaticity breaks down and the fluctu- 
ations are frozen at a finite value. The inverse of the 
nonlinear coupling defines a minimum time = 1/kN 
that is needed to obtain squeezing. Only for extremely 
long times the fluctuations approach An = 0. 

Thus high number squeezing requires very long split- 
ting times, and a linear improvement in squeezing is 
achieved only by an exponentially longer splitting time. 
For typical parameters, the necessary splitting time 
might even exceed the condensates coherence and life 
time. 



hnear interaction twists the distribution around the z- 
axis, associated with the action of (one-axis twisting 
model [36 ). The twist rate is zero at the equator, in- 
creases linearly when moving towards the poles, and the 
twist direction is opposite on the north and south hemi- 
sphere. For an initial binomial state, the nonlinear cou- 
pling twists the state such that it becomes a squeezed 
state, but along a direction which is not the equator of 
the Bloch sphere. In general, the number squeezing is 
low in this case. However, when a constant tunnel cou- 
pling is added [see flrst term in Eq. ([5])], it has the effect 
of rotating the spin squeezed state around the x-axis. By 
varying the strength and time of the tunnel coupling, an 
optimal number squeezing can be achieved. These opti- 
mal parameters are given by analytical formulas [35] . 

For applying this to the splitting process, we use a 
control which consists of two stages: 

^^(t) = Oo(^l-^^e-*/*^+Q,. (7) 

First, we exponentially turn off Q{t) with a time constant 
tc, and in the second stage the control is kept constant. 
Thus, the longer tc the more the state is adiabatically 
squeezed in the first stage. It turned out that a more 
adiabatic splitting allows for more squeezing in the end, 
but at the price of a longer squeezing time. Optimiz- 
ing for these two parameters for given time intervals, we 
find much better squeezing than with a pure exponential 
control. Results for different are displayed in Fig. |3] 



D. Optimal control theory 

Optimal control theory (OCT) ^E2lH provides an 
alternative for achieving efficient number squeezing on 
short time scales. OCT is a mathematical device which 
allows the optimization of a control objective (such as 
high number squeezing) for a system whose time evolu- 
tion is subject to a constraint (e.g., time evolution gov- 
erned by the two-mode model). In general, the "opti- 
mal control" is obtained numerically through an iterative 
procedure [3X 44 • We start by defining a cost function 
which quantifies the optimization target. For the case of 
number squeezing it reads^ 

j(c, Qs) = (c(r)i j,2|c(r)> + J r{ns)^dt , (8) 



C. Two-parameter optimization 

Law et al. [34] suggested an alternative control scheme 
for achieving a much higher number squeezing in com- 
parison to the exponential switching-off. This is achieved 
by applying for a given time interval a constant tunnel 
coupfing. To understand the working principle of this 
control, we note that in the Bloch-sphere picture the 
tunnel coupling corresponds to a rotation around the x- 
axis, associated with the action of Jx, whereas the non- 



where T is the terminal time of the control sequence. 
The second term, which is weighted by the regularization 
parameter 7, penalizes rapidly varying tunnel coupfings 
and is needed to make the control problem well posed. 
We choose the value 7 = 10~^ throughout this work, such 



^ In H we additionally make the replacement Q 
ensure positive tunnel couplings. 



■ fi? in order to 
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that the first term dominates. J(C, Qs) is now optimized 
under the constraint that C obeys the time evolution of 
the Schrodinger equation, with the Hamiltonian given in 
Eq. ([5]). This is achieved by introducing a Lagrange 
multipher C for the constraints. We define a Lagrange 
function 



L(C, C, Qs) = J{C, Qs) + Re( C, iC - HC 



(9) 



where C is also called the adjoint variable. The bracket 
stands for (u,v> = {u{t)\v{t)) dt. ^(0,0,0^) has 
a saddle point at the minimum of J{C,Qs)- Thus, the 
Frechet derivatives with respect to C, C and Qs must 
all be zero for the optimal control. From this we obtain 
Euler-Lagrange equations, together with the initial and 
terminal conditions, 

2|C> = H\C), C(0) = Co (10a) 
i\C) = iJ|C>, i\C{T)) = 2J^\C{T)) (10b) 
= 205Re(C|4|C> . (10c) 

The equation in the third line is subject to the (desired) 
boundary conditions 0^(0) = Qg and ^^(r) = and it 
is fulfilled for the optimal control. For any other, non- 
optimal control the third equation does not hold. In- 
stead, the so-called control equation 



VJ = + 20,Re(C|4|C) . 



(11) 



applies, which can be interpreted as the gradient of the 
cost functional with respect to a variation of the control 
field. It allows for the implementation of an iterative 
numerical scheme to determine the optimal control. 

One can thus use the following algorithm to iteratively 
improve the control: For an initial guess, one first cal- 
culates the forward equatio n Eg . ( |10a| ), and thereafter 
the backward equation Eq. (10b). From the gradient in 



Eq. (11), a new search direction is found, from which 
an optimization scheme, such as the nonlinear conju- 
gate gradient or quasi-Newton one [HI [45], determines 
a new control. For the time propagation we use a Crank- 
Nicolson scheme similar to Ref. [T9^. 

The direct implementation of the above optimization 
scheme can lead to problems. One of them is the appear- 
ance of discontinuities of 0(t) at the initial and final time. 
Another one is that the algorithm is sometimes unable 
to optimize the global structure of the control but sticks 
very close to the initial guess. Both problems can be 
circumvented by using a more refined approach gH]: in- 
stead of the usual L^-formulation of the control problem, 
we employ a control XeH^, i.e., in the space of func- 
tions with weak derivative again in L^. The Lagrange 
function (and all inner products appearing in the opti- 
mization algorithm) has then to be reformulated using 
the inner product {u,v)h^ = {u,v)l'2. In particular, 
we get ^{X,X)h^ for the regularization term. While the 
forward and adjoint equations remain unchanged within 



TABLE I: (Color online) Definition of Josephson and Fock 
regime used in the text. There exists also the Rabi regime 
{Q/k ^ N), which, however, does not play a role in the 
context of number squeezing. 



Regime 


Criterion 


Characteristics 


Josephson 




Competition between 






Q and k 


Fock 


n/K < i/N 


Non-linear interaction 






K dominates 



this approach, the gradient is given now by a Poisson 
equation 

- ^VJ = -^Qs - 2Q,Re(C|4|C> , (12) 

which distributes local changes of the rhs to the com- 
plete time interval. Since second time derivatives of Qs 
appear on both sides of the equation, the gradient has 
the same degree of smoothness as the control. Further, 
the final and terminal conditions of the control are ful- 
filled exactly. In the context of the GP-equation it was 
demonstrated that this approach leads to a smoother 
control and quicker convergence than the common one 

m- 

While the standard approach is sufficient for our op- 
timizations within the generic model, we will find clear- 
cut advantages of the formulation in case of the real- 
istic simulations. 



E. OCT Results 

In a bosonic Josephson junction one distinguishes be- 
tween the Rabi-, Fock- and Josephson regimes [20[ l39]. 
which are characterized by different values of Q/k listed 
in Table [T] In the Josephson regime there is a competi- 
tion between tunneling and the non-linear interaction k, 
while in the Fock regime k dominates. It turns out that 
our OCT fields decisively depend on the initial guess for 
the tunnel coupling 0(t), although the final states have 
similar number squeezing and are both suited for atom 
interferometry purposes, as will be discussed in section |Vl 
We essentially find two types of solutions, which can be 
classified according to physical regimes discussed above 
as Josephson and Fock solutions. 

For the Josephson solution, the initial guess (dashed 
lines in Fig. [4| is chosen such that after a first exponen- 
tial decrease of Q the control decreases linearly. In this 
second stage, the system is in the Josephson regime with 
a competition between tunnel coupling and nonlinear in- 
teraction. With this initial guess, OCT comes up with 
an optimized Q that shows oscillations in the Josephson 
regime, as shown in Fig. [4| Also the fiuctuations An os- 
cillate, and at the final time the system is brought to a 
highly squeezed state. 
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FIG. 4: (Color online) Typical results from Josephson OCT 
for the generic two-mode model for different splitting times 
and (a) = 100, (b) = 1000. The initial (optimized) 
Q is shown by the dashed (solid) line in the lower panels of 
the figures. The upper panels of the figures show the number 
fluctuations. The bright lines show estimates from a para- 
metric oscillator model, see Sec. IIIF (a) The Bloch spheres 



demonstrate the state for different times. 



For the Fock solution the initial guess is a simple ex- 
ponential decay (dashed lines in Fig. [5|. Here, the op- 
timized displayed in Fig. [s] is most of the time very 
small (and thus can be classified to be of Fock type). 
Only for two short time periods it rises to larger val- 
ues. Comparison between exponential, Josephson, and 



Fock OCT sphtting is made in Fig. [T2j Best squeezing is 
achieved within the last approach. 



F. Interpretation 

Josephson optimization. — We next interpret the oscil- 
lating behavior in the Josephson control of Fig. |4| To 
this end, we exploit that for large atom numbers the 
state vector C can be approximately treated as a contin- 
uous variable [19]. We consider C{k) as the atom- number 
wave function and k as the number difference between the 
left and right well. The matrix elements of Eq. Q are 
Taylor expanded in powers of and C{k ± 1) is Tay- 
lor expanded around k. We then obtain a Schrodinger 
equation for the variable C{k) 



iC{k) = 



n{t)N 
^dk^ 



+ 



n{t) 



+ 2«: /c 



C(/c), (13) 



which can be identified as a parametric oscillator in the 
variable k [47j . One could think of the parametric oscilla- 
tor as a harmonic oscillator with constant frequency but 



FIG. 5: (Color onhne) Typical results from Fock OCT for 
the generic two-mode model for (a) TV = 100, (b) = 1000. 
The initial (optimized) Q is shown by the dashed (solid) line 
in the lower panels of the figures. The upper panels of the 
figures show the number fluctuations, (a) The Bloch spheres 
demonstrate the state for different times. 



with a driving term proportional to k"^. Thus, paramet- 
ric amphfication translates to a nonhnear driving process, 
where the fluctuations (and not the mean atom number 
difference) are subject to parametric amplification. 

In appendix|A|it is shown that if the oscillator is driven 
by an oscillating tunnel couphng Q{t) = + cos (cut) 
(we assume <C Qq) with approximately twice the res- 
onance frequency tUres = f^N Oq, the envelope for the 
width of the initial number distribution C{k) decays ex- 
ponentially. For the envelope of the number fluctuations 
we obtain |i48] 



(An) envelope ~ AnQ CXp 



2(0o + kN) _ 



(14) 



The bright hues in Fig. |4| show results of this model: 
the oscillation period and shape of the control indeed 
mimic ^ocrii), and the envelope of the number fluctu- 
ations is in good quahtative agreement with AnocT(^)- 
The fluctuations are reduced to values of approximately 
An/Ano ~ 0.3 — 0.5, depending on the number of atoms. 
During the flnal stage of the control, nonlinear cou- 
pling dominates over tunneling, and the approximations 
introduced above fail. Squeezing in this interaction- 
dominated regime can be understood in terms of the one- 
axis twisting model discussed above. 

Fock optimization. — The physics underlying Fock opti- 
mization can be regarded as a generalization of the model 
underlying the two-parameter optimization [34] . The re- 
sults from Fig.[5|can then be interpreted as follows. After 
the exponential splitting, we end up in a squeezed state 
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which, however, is not a number-squeezed state. The first 
O pulse after the exponential decay rotates the squeezed 
state towards the equator, but a little tilt remains, which 
leads in the ensuing time evolution to squeezing caused 
by the non-linear interactions. During this stage Q is al- 
most zero. In the end, a second O pulse rotates the state 
to the equator, such that the number squeezing becomes 
maximized. 

In Fig. [3] we compare the two-parameter optimiza- 
tion with Fock OCT. In a wide regime of splitting times 
and atom numbers, the squeezing achieved within both 
schemes agrees, except for small atom numbers and long 
splitting times, and large atom numbers and short spht- 
ting times. This is due to the fact that for the constant 
control scheme [34^ there exists an optimal squeezing and 
an optimal squeezing time, which are both directly pro- 
portional to [35J. In case of OCT however, we are not 
restricted to a single control strategy, such as the mere 
value. As a consequence, OCT can come up with more 
improved control strategies and performs better than the 
more simple schemes. 



IV. OPTIMIZATION OF A REALISTIC MODEL 

We next proceed to a more reahstic model for the spht- 
ting process, which captures also the spatial dynamics 
and relies on the control parameter used in experiments, 
related to the double well separation of the trapping po- 
tential. As consequence, we can also control the oscil- 
lations of the condensate when it is split fast. Methods 
for including the orbitals describing the spatial dynamics 
rely on variational principles, as for example in Ref. [22] 
where a Gaussian ansatz for the orbitals was used. 

In our approach we employ the multi configurational 
time dependent Hartree for Bosons method (MCTDHB) 
[33l l49] , where the orbitals are included self-consistently, 
without any assumptions about their shape. Although 
the MCTDHB framework is formulated for an arbitrary 
and fixed number of orbitals, we will stay with two, as 
they are sufficient to account for condensate splitting. 



A. Doublewell potential 

In our calculations we consider condensates in elon- 
gated magnetic traps with a tight transverse and a weak 
longitudinal confinement as can be manipulated easily 
on atom chips. For the transverse confinement, we use 
the 2-d potential given in Ref. [14]. In the longitu- 
dinal direction we take an harmonic confinement with 
trap frequency (jJ\\. Condensate splitting is achieved by 
modifying the trap potential along a transverse direc- 
tion, through the confinement potential discussed in Sec. 
II B and Ref. [14]. The trap is then split along the 
x-direction by varying the radio- frequency field B^i{\), 
where A = — | corresponds to the unsplit trap which can 
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FIG. 6: (Color online) Inverse adiabatic squeezing time scale 
UqN versus ID-interaction strength Uq for several A^, where 
each point along a straight line corresponds to a different 
aspect ratio (connected by black lines to points with same 
aspect ratio), as given in the figure. The straight lines are 
guides to the eye. 

be approximated by an harmonic well with a transverse 
trap frequency uj±_ = 27r • 2 kHz. 

In our model we assume that the splitting direction is 
decoupled from the other ones. For typical aspect ra- 
tios the dynamics in the weakly confined lon- 
gitudinal direction is much slower than in the splitting 
direction. The 3D interaction strength is given in a con- 
tact potential approximation Uq^ = Anh'^a/m, where a 
is the s-wave scattering length and m the mass of an 
atom. For calculating the ID interaction parameter Uq^ 
we consider the stationary CP equation ^SSj in 3D with 
the confinement potential v{r) = Vx{x)-\-Vy{y)-\-Vz{z) and 
define the single particle Hamiltonians hi = —^W'^^Vi 
{i = x,y,z). We assume that the wave function fac- 
torizes, 0(r) = (f)x{x)(l)y{y)(f)z{z), and we integrate out 
the y and z coordinates, (px then satisfies a non- linear 
Schrodinger equation with the Hamiltonian 

Hx = hxc^x + U^'^iN - l)WyWz\c^x\^c^x , (15) 

where Wy = J dy\(j)y{y)\'^, with an analogous expression 
for Wz- 4>y and (pz satisfy similar equations with Hy and 
Hz, respectively. Employing imaginary time-propagation 
simultaneously for (/)x) (t^y and 0^, we obtain the ground 
state of the trap. The ID interaction parameter for spht- 
ting in x-direction is then given as Uq = UQ^WyWz- 

In Fig. [6] we plot UqN, which determines the squeezing 
time scale for the quasi-adiabatic exponential splitting, 
versus the 1-d interaction parameter Uq. The straight 
lines correspond to a fixed atom number taken for 
various values of the aspect ratio of the trap uj±/uj\\ = 
50, 100, 300, 1000. For small aspect ratios Uq is largest. 
It is apparent from the figure that neither Uq nor UqN do 
increase linearly with the atom number A^. Instead, for 
higher the wave function spreads out more along the 
longitudinal direction, and thus Uq is reduced. Therefore, 
the value UqN = 1 SiS assumed in most of our optimiza- 
tions is quite realistic. The chemical potential is between 
/i ~ 27r • 2 kHz (large aspect ratio) and /x ~ 27r • 4 kHz 
(small aspect ratio). 
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B. MCTDHB 



1. Equations 




The MCTDHB equations are obtained by making an 
ansatz for the many-body wave function in a basis of two 
orbit als, which are ahowed to depend on time, and apply 
variational calculus. In the following we summarize the 
MCTDHB equations [33] . Since our trapping potential is 
symmetric, the orbitals exhibit gerade {(/)g) and ungerade 
(0e) symmetry. 

The orbitals equations read 



;o 10^ 



(b) 



i(t)g 



= P 
= P 



,(16) 



/l0e + (/f|0,P+/||0eP)0e + /< 




The coefficients are given by 

fk = Uo{p}^^pkkkk , fk = ^Uq{p}^1 pkqkq , 
fk = Uo{p}^^pkkqq , 



(17) 



where k is either g oi e, and q the opposite. The one- 
and two-body reduced densities [50] pkq and pkqim, re- 
spectively, are given by 



X(G) 



FIG. 7: (Color online) (a) Mean population difference a — 
{pgg — pee)/N of the orbitals normahzed by (a is the co- 
herence factor) and number fluctuations , (b) nonlinear coef- 
ficients of the orbitals equations from Eq. ( |17| , (c) the tunnel 
coupling, and (d) the two-body matrix elements of Eq. \2l\ . 
Parameters are UqN — 1 and = 100 (in the following fig- 
ures we always use these values, unless stated otherwise). 



Pkk = {C\a\ak\C) , pkkkk = {C\alakalak\C) , (18) are the two-body matrix elements. Eqs. ([T6| and ([To]) are 



Pkqkq = {C\alakalaq\C) , pkkqq = {C\alalaqaq\C) . 

Note that because the orbitals have different parity, the 
one-body reduced density is diagonal and all matrix el- 
ements of the two-body reduced densities with uneven 
combinations of indices vanish [42 . We distinguish be- 
tween two types of nonlinearities: the terms of type 
and are GP-hke, where the orbitals experience 
a real potential proportional to the moduli of the or- 
bitals. The terms of type fk are complex contributions, 
and couple the two orbitals in a more complicated fash- 
ion. P = 1 — \(pg){(l)g\ — |0e)(0e| is a projector which 
guarantees that the orbitals stay normahzed throughout. 
The number distribution vector obeys the equation 



now coupled equations which have to be solved simulta- 
neously. 



2. Ground state properties 

We obtain the ground state of the trap for a fixed 
splitting A by applying imaginary time-propagation to 



Eq. (16) and simultaneously diagonahsing Eq. (19). This 



ot 



The two-mode Hamiltonian is given as 

1 
2 



(19) 



(20) 



k,q,l,m 



where the prime indicates that the sum only runs over 
even combinations of indices, which is due to the par- 
ity of the orbitals. We have introduced annihilation 



(t) 



and creation operators 
ade orbitals. Further, we have Jx = ^(a^ac, 

n = {(l)Sh\(l)e) - {(l)g\h\(l)g 



for the gerade and unger- 
2 g'^g cigO^e) and 
The integrals 



Wkqim = Uo j dx(t)l(x)(t)q(x) 



(21) 



gives us the MCHB groundstates. 

Consider first the unspht trap: The gerade orbital is 
populated to more than 99% [Fi g.[7] (a)]. The only signif- 
icant nonhnear couphng in Eq. ( [l6{ is due to f^^UoN 
[Fig. [t] (b)]. This reflects that the ground state of a 
harmonic trap is well described by the CP equation, 
and the quantum depletion out of the condensate is low 
[42] , At the same time, the ungerade orbital is affected 
strongly, through the coefficient /g, such that the energy 
gap between the orbitals becomes larger in comparison 
to the harmonic oscillator states. For increasing interac- 
tion strength the gap decreases, and the excitations to 
the ungerade orbital are enhanced. We note that the 
structure of the term proportional to fe is reminiscent 
of the Bogohubov-de Gennes equations, which describe 
weak excitations out of the condensate [39j . 

When the potential is split into two separate wells, the 
energy splitting decreases towards zero, and so does the 
tunnel coupling [Fig. [t] (c)]. In addition, the ungerade 
orbital becomes more and more populated, as shown by 
the red dashed line in Fig. [t] (a). For the split trap, 
the wave function ^^^g = (0l ± (/)r)/V^ can be decom- 
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posed into left and right orbitals, which are solutions 
of Gross-Pitaevskii equations with the nonhnearity for 
~ N/2 atoms. In Fi g, [t] (d) we plot the two-body matrix 
elements from Eq. ( |21| ). They drop by a value of about 
two during sphtting, and are all equal only for the spht 
trap.-^ 

Number fluctuations in the spht condensate are shown 
in Fig. [t] (a). The observable for detecting an atom on 
the left side (i.e., in the left well of the split trap) reads 



Nl = dx^"^^ = N/2 + dd)gae + (falcLg , 

J —oo 



(22) 



where d = (f)*g{x)(f)e{x) dx is the overlap integral over 
the negative half axis. A similar expression is obtained 
for Nr by restricting the integration to the positive half 
axis. The number difference operator is then given by 
AN/2 = {Nl - Nr)/2 = ddlde + d*dldg, where we have 
exploited the symmetry of the orbitals. For the split trap, 
where (f)g^e = (0l ± 4'r)/V2 holds, we obtain d = ^ and 

AN/2 corresponds to Jz = (a|^aL — «i?«i?)/2, otherwise 
d < ^. For a dynamic splitting of the trap, d can acquire 
an additional phase, which, however, is irrelevant when 
computing any expectation value. 



C. Optimality system 

The cost functional consists of a squeezing and a reg- 
ularization term, analogous to the cost function for the 
generic model in Eq. ([sf . Additionally we require that the 
spatial dynamics is stationary at the final time to avoid 
oscillations of the condensates in the split trap [32]. It 
turns out that this can be achieved by trapping the or- 
bitals in the GP ground and first excited state, respec- 
tively, which become the "desired" orbitals 0^ and 0f 
of our optimization approach. We find that, at least for 
moderate interactions, these orbitals are stationary, inde- 
pendent on the precise value of the number fiuctuations. 
We then have 



j{4>„<t>e,c,x) = |L(c(r)|^|c(r)) 



+?[2-|(0s(r)l<>r-|(0e(r)|0^>| 



-I 



[m]^dt, 



where 71, 72 and 7 weight the different control objectives. 
For the control A we choose fixed initial and terminal 
conditions, i.e., we require an initial unsplit trap and a 
final split trap. The optimality system for MCTDHB, 



which is quite involved, is derived analogous to Sec. |IIID 
Details are given in Appendix [B] 



D. Numerical implementation 

The MCTDHB equations are highly nonlinear equa- 
tions, due to the orbitals nonlinearities, the state vectors 
dependence on tunnel coupling and two-body matrix el- 
ements, and the projectors. For the time integration we 
use a Modified Crank-Nicolson method with Newton 
iterations at each time step. Details are given in ap- 
pendix |C] Using an implicit scheme is very advantageous 
in connection with OCT, since preservation of the norm 
is crucial. This is achieved here with finite step size and 
accuracy 0{At^). For the backward equations, we can 
use the wave function from forward integration due to 
the constant time step. 

A very useful criterion for testing the correct imple- 
mentation of the optimality system, which is usually 
quite involved, is as follows. We choose test functions 
u and V, and compare the finite difference {J(u + ev) — 
J{u — £v))/{2£) to the inner product (v, VJ(n)). If it 
converges for sufficiently small e, the gradient is indeed 
a viable search direction. If not, there is an error in the 
implementation of the optimality system. At the same 
time, this procedure is a test for the accuracy of the nu- 
merical scheme. 

For the OCT results within MCTDHB, the function 



space discussed in Sec. HID turned out to be of impor- 



(23) 



tance. Optimization within yields controls related to 
Josephson optimization, whereas within we recover 
Fock optimization. Within predominantly the local 
structure of the control is modified, and oscillations are 
added. These also lead to condensate oscillations, which 
are hard to stop at the end of the control sequence and 
thus confiict with the trapping term in the cost function. 
Therefore, the choices of the best initial control fields and 
the ratio 71/72 turned out to be crucial. In contrast, the 
optimization is capable of changing the global struc- 
ture of the control, and is therefore less susceptible to 
such details. 



E. Limitations of a two-parameter optimization 



In our comparison of the generic two- mode model and MCHB, 
we neglect the effect of different Wkqim before and after splitting, 
which would lead to a slight renormalization of the tunnel cou- 
pling ^11^. We use the value of Wgggg of the split trap, which 
deviates by a factor of approximately two from the value for the 
unsplit trap. Considering that the two-body matrix elements of 
left-right and gerade-ungerade orbitals differ by a factor of two, 
we are left with k = Uo/2. 



In principle, a two-parameter optimization, as pre- 
sented in section [TTT| is also possible within the realistic 
model. However, the splitting time scale tc cannot be 
made arbitrarily small due to the spatial orbital dynam- 
ics. For a too fast splitting, the condensates and thus also 
the tunnel couphng strongly oscillate. Although we found 
that small oscillations in Q, which are always present, do 
not disturb much, strong oscillations typically prohibit 



10 




has been achieved, say to a value of A = 1.2, which corre- 
sponds to a double well separation of about 3 fim. Dur- 
ing this separation the number fluctuations might still 
change, which makes a simple control strategy like the 



Time 



FIG. 8: (Color online) (a) Trapping cost function for a two- 
parameter control for several UqN and tc- (b) Example for 
the density and the tunnel coupling for a specific example 
{UqN — 1, tc — 0.23, marked by the red point), which both 
strongly oscillate. Here we have = 100, but results are 
similar for difi^erent A^. 



large squeezing. In addition, too strong condensate os- 
cillations during splitting lead to unwanted excitations 
and to the breakdown of the two-mode approximation. 

In order to quantify those oscillations we calculate the 
cost functional for trapping from Eq. ( [23| ) for different tc 
and Uq (the value of does not matter for the present 
concerns), see Fig. [sja). In panel (b) we report the 
atomic density and the tunnel couphng for a specific ex- 
ample. At these parameters the cost functional is larger 
than say 0.025, O strongly oscillates, and no significant 
squeezing can be achieved. 

From this we conclude that the two-parameter opti- 
mization works only for large enough splitting times, 
which is restrictive primarily for stronger interactions 
and larger particle numbers. The green lines in Fig. [3] 
report the shortest times on which the two-parameter 
control works within MCTDHB, where the shortest al- 
lowed time scale tc for each UqN is taken from Fig. [8| 
and two-parameter optimization performed within the 
generic model to get an estimate. For easier comparison, 
all results in Fig. [3] are rescaled such that the adiabatic 
squeezing time scale {I/UqN) is the same. 

Decoupling the condensates is necessary after splitting 
to avoid any residual tunnel couphng. Thus we sepa- 
rate the split condensates further after high squeezing 



two-parameter optimization from Sec. Ill rather prob- 
lematic. This can be easily accounted for within OCT as 
we discuss below. 



F. Optimization and results 

Josephson optimization. The initial guess is obtained 
from the tunnel coupling for the ground states shown 
in Fig. [tJc). We choose the ramp such that the tunnel 
coupling for the static orbitals first quickly drops in an 
exponential fashion, and add a linear function, similar to 
our optimization within the generic model. In the first 
stage, the trap is split quickly up to the point where the 
tunnel couphng is comparable to the nonhnear interac- 
tions. In the second stage, the interplay between tunnel- 
ing and nonlinear coupling can be exploited to bring the 
fluctuations down. 

Our optimization is performed sequentially, where (i) 
we optimize squeezing and bring the orbitals to a halt 
(here both 71 and 72 are finite), and (ii) further separate 
the condensates and trap them (here only 72 is finite). 
After these two steps (hi) a waiting period follows (no 
optimization here), which mimics the stage where the 
condensate in one well of the interferometer might ac- 
quire a phase. With this initial guess, OCT again comes 
up with a control reminiscent of parametric amplifica- 
tion, as shown in Fig. [9| The large oscillations are due 
to the parametric oscillations, and the small ones due to 

Iwe 



the unavoidable oscillations of the orbitals. In Fig. 10 



show the tunnel coupling of the control, which exhibits 
the oscillatory behavior expected for parametric amplifi- 
cation [30 . When comparing to the generic model, one 
expects nonlinear renormalizations to O due to the over- 
lap integrals Wkqim [51^, which heavily oscillate here. In 
this sense, the comparison of Vt in both models can be 
only a qualitative one and a negative has no physical 
significance. 

The parameters for which Josephson optimization 
works are restricted by the time scale of the squeezing 
dynamics, determined by 1/UqN ^ which should be much 
larger than the time scale determined by the transverse 
trap frequency. For UqN ^ 1 the strong oscillations of 
the orbitals would spoil the control. Stability due to un- 
certainties in the pertinent parameters is also displayed 
in Fig. [9] The control is stable to about 10% variations 
of the atom number or the nonlinearity. Changes in trap 
frequency or shape of the double-well, on the other hand, 
infiuence the oscillations in the orbitals more significantly 
and can degrade the control, as discussed in the figure 
caption. 

Fock optimization. — When using the optimization 
the control becomes smoother. The tunnel couphng de- 



picted in Fig. 10 shows a behaviour characteristic for 
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FIG. 9: (Color online) In the left panel, a typical solution 
for Josephson optimization for = 100. Changing interac- 
tions (red line) or atom number (green line) by 10 % do not 
appreciably change results. However, the control scheme is 
sensitive to the (transverse) trap parameters. Here, the de- 
viation of the initial harmonic trapping frequency should not 
exceed 0.1% (dashed blue hue) and is crucial for 1% (sohd 
blue line). In the right panels, the corresponding orbitals are 
shown. 




FIG. 10: (Color onhne) Tunnel couphng within MCTDHB 
for the Josephson (dark blue line) and the Fock (bright red 
line) solution from Fig. [9] and Fig. 11 (a), respectively, where 



the second one is scaled by a factor of 10. 



Fock optimization. Typical solution are given in Fig. [TT 
In contrast to Josephson optimization, suitable values of 
7i and 72 can be found more easily, and the condensates 
always end up in a stationary state. Here, we demon- 
strate that we can separate the condensates quite far, 
while they are at rest at the end of the control sequence. 
When separating less, the condensates motion during the 
control sequence is even less oscillatory [53]. The sensi- 
tivity to the (transverse) trapping potential is much less 
severe than for Josephson optimization. 

Within Fock OCT we find solutions also for stronger 
interactions {UoN = 5), where the orbitals dynamics is 
more crucial, as shown in Fig. 11 (b). 



Squeezing values for both OCT strategies in the regime 
of weak interactions {UqN = 1) are shown in Fig. [12] 
where we also take into account the time it takes to un- 
couple and trap the condensates. We also present a com- 
parison with the generic model. 



V. OPTIMAL SQUEEZED STATES IN ATOM 
INTERFEROMETRY 

In atom interferometry [3], a phase difference between 
the split condensates is acquired due to the interaction 
with a weak external potential. For the readout of the 
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FIG. 11: (Color online) In the left panels typical OCT solu- 
tion for Fock optimization for N — 100 with (a) weak interac- 
tions UqN — 1 and (b) strong interaction UqN — 5. Exponen- 
tial (OCT) control and squeezing is displayed by the dashed 
(solid) lines. In the right panels, the orbitals are shown. 
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FIG. 12: (Color online) Comparison of squeezing (symbols) 
within the generic model and MCTDHB for UqN — 1, com- 
pared to quasi-adiabatic exponential splitting (gray lines). 
The latter is calculated within the generic model, since MCT- 
DHB results are very similar. Best results are obtained for 
Fock optimization and large TV. 



phase difference, several schemes exist. In time-of- flight 
(TOF) experiments the trap is switched off, the conden- 
sates expand freely and finally overlap. From the in- 
terference fringes the phase can be deduced O [HI fTT] , 
The densities for typical control sequences consisting of 
sphtting, waiting (or phase accumulation) time, and re- 
combination in TOF, are shown in the upper panels of 

Fig.g 

Another method is phase sensitive recombination of 
the split condensates, as has been demonstrated in [54] , 
Even more desirable would be a Mach-Zehnder (MZ) like 
setup [20 , where a 7r/2 pulse of the tunnel coupling [55] 
transforms the initial number squeezed state into a phase 
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FIG. 13: (Color online) In the upper panels, the density for 
(a) an exponential and (b) an optimized splitting is shown. In 
the last sequence the trap is switched off and the condensates 
start to expand and overlap (TOF). In the lower panels, the 
interference in momentum space, obtained from the Wigner 
function W{x — 0,p) [53 , is plotted. When the condensates 
are released, the interference in momentum space is trans- 
formed to an interference in position space. The fringe con- 
trast is reduced in case of exponential splitting due to phase 
diffusion, and squeezed states are much more robust. For a 
perfect fringe contrast, the limitations to the phase sensitivity 
are given by the quantum noise. 



squeezed state, and after phase accumulation the phase 
information is transformed into number information by 
another 7r/2-pulse. The readout is achieved then by a 
relative atom number measurement [18 . 

Squeezed states have reduced quantum fluctuations 
and thus phase (number) squeezed states have the poten- 
tial to increase the phase sensitivity of an atom interfer- 
ometer below the standard quantum limit Ac/) = l/\fN 
[H [TB] in a TOF (MZ) setup. Phase squeezed states can 
be obtained from number squeezed states by applying 
a 7r/2 tunneling pulse. The phase sensitivity in terms 
of shot noise (I/a/TV) is given by the factor of useful 
squeezing = [20l [21], where a is the coher- 

ence factor [56]. It determines the fringe contrast seen 
in TOF experiments and is given by the polarization 
of the state along the x-direction of the Bloch sphere 



a = 



— {pgg~Pee)/N . In principle, its value tends to 
zero for very small number fluctuations due to the uncer- 
tainty relation between number and phase [see Fig. 7 (a)], 
but for the squeezed states of our concern here, the hm- 
iting factor is the phase diffusion to be discussed below. 

is bounde d from below by the fundamental Heisenberg 
hmit U = ./2/N PL. 




Time 



FIG. 14: (Color onhne) for optimized (symbols) and expo- 
nential splitting (gray lines). The quasi-adiabatic exponential 
splitting is shown for several time scales for TV = 100. After it 
reaches a minimum, it rises again due to loss of phase coher- 
ence. With optimized splitting, lower is obtained at shorter 
time scales. The symbols are the same as in Fig. The 

results are compared to the Heisenberg limit — \/ (hor- 
izontal hues). For larger atom number = 500 (MCTDHB) 
and = 1000 (generic model) optimized splitting yields lower 
^R, whereas for adiabatic splitting the curves are similar. 
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FIG. 15: (Color online) Coherence factor a for typical 
Josephson OCT (dark red) and Fock OCT (bright green) so- 
lutions, where it decays slower than for the quasi-adiabatic 
exponential splitting (broken line). The Josephson solution 
shows up oscillations in the phase coherence during splitting, 
because the driving of the parametric oscillator affects both 
conjugate observables of number and phase. 



In Fig. [T4| we plot (r for the optimal and the quasi- 
adiabatic exponential splitting, demonstrating the im- 
provement gained by applying OCT. Note that for the 
exponential control the minimum of $^r is reached at fi- 
nite values of O, where the squeezing still changes in time. 
Finite tunnel couphng Q leads to phase locking between 
the condensates [T71 157] , which comphcates measurement 
protocols. Thus, additional control strategies would be 
necessary to uncouple the condensates at the right time. 
Josephson and Fock OCT results yield quite similar ^r, 
since in the first case there is less squeezing but better 
phase coherence. 

In presence of interactions, the condensates are sub- 
ject to phase diffusion [4 which quickly reduces the phase 
sensitivity. On the Bloch sphere, the states are twisted 
around the z-axis, which leads to a decrease of the phase 
coherence. Whereas phase diffusion, which is propor- 
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tional to An 0], is very rapid for binoraial or phase 
squeezed states, number squeezed states are much more 
robust. However, using number squeezed states in a TOF 
interferometer comes at the price of a phase sensitivity 
above shot noise. In Fig. [T5| we compare the phase co- 
herence for various splitting protocols during and after 
sphtting. Obviously, both OCT results are much more 
stable against phase diffusion as compared to an expo- 
nential splitting on the same time scale. This leads to an 
improved fringe contrast, as is shown in the lower pan- 



els of Fig. [T3| using the Wigner function. Thus, squeezed 
states allow for much longer phase accumulation times. 



Conclusions and Summary 

We discussed optimal control of the splitting of a Bose- 
Einstein condensate such that number squeezing is mini- 
mized. For a physically realistic modehng of the splitting, 
we employed the MCTDHB equations where the spatial 
dynamics is included and the control is through the trap- 
ping potential. This allows for a direct application of our 
sphtting protocols to experiments. Optimal control of 
MCTDHB turned out to be sensitive to different math- 
ematical strategies, which give rise to different physical 
solutions. The optimization of a generic two-mode model 
allowed us to obtain an intuitive understanding of the un- 
derlying physical processes. Two control strategies have 
been identified, which were compared to a more simple 
two-parameter optimization. The latter, however, has 
only limited validity in case of realistic splitting. Re- 
sults have been given for several squeezing time scales 
and interactions strength, demonstrating that with re- 
fined protocols number squeezing is achieved at least one 
order of magnitude faster than with exponential split- 
ting. Finally, we analyzed the squeezed states in context 
of atom interferometry. 
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^ The question which atom number is ideal for reducing phase 
diffusion in the traps under consideration can be answered as 
follows. The phase diffusion rate |19J is determined by ~ UqAti. 
Prom Fig. [g] we infer that larger N means less interactions Uq, 
and in addition the relative squeezing achieved within OCT is 
better for larger N (in contrast to the absolute squeezing An). 
Finally, it turns out that the phase diffusion rate is roughly the 
same for all N (when using OCT splitting). However, higher 
number squeezing means enhanced phase fluctuations. 



APPENDIX A: DERIVATION OF THE 
PARAMETRIC OSCILLATOR MODEL 

In this appendix we present details on how to derive 
the resonance frequency for parametric amplification and 
Eq. p4| ). We start by considering the Hamiltonian of 
Eq. ( |13[ ), which for a? = is a harmonic oscillator with 
a time dependent mass m{t) = q^jj^ and frequency 
uj{t) = Q{t). For the further steps it helps to rewrite 
the Hamiltonian using normal variables, which read [47] 



1 



We obtain 



iC = 0(t){(ata + -) + ^(a^ + a''^)}C , 



(Al) 



(A2) 



with Cl{t) = Q{t) -\- kN. In order to get rid of the time 
dependence in the noninteracting part, we transform to 
the variable 



sit) = / dt'n{t') 



(A3) 



We are now left with an harmonic oscillator which is 
driven nonlinearly, a model which is frequently used, e.g., 
in quantum optics [48]. In an interaction representation 
with respect to Hq = a + ^ the Hamiltonian reads 



H = 



2Q{t{s)) 



(a^e-^- 



(at)2e2-) 



(A4) 



In the spirit of our OCT control (see Fig. |4|) we assume 
that the tunnel coupling is an oscillating function with 
some constant (or slowly varying) offset. The resonance 
condition is fulfilled for ^}{t{s)) = + cos 2s, where 
the dependence of t on s for this case can be obtained 
from Eq. (IaI: 



: arctan 



(1 - /3)tans 



■}. (A5) 



Assuming = 75-^ <C 1, we have to a good approximation 
simply s = Clgt. This gives us the resonance frequency 

For discussing the time scale of the decrease of the 
number fluctuations, we separate the Hamiltonian of 



Eq. (A4) into a resonant and a non- resonant part. Us- 
i^g i+/3cos\2ujrest) ^ 1 " COS (2u;^es0: the resonant part 
reads in a rotating wave approximation [48] 



(A6) 



Solving for the fluctuations of a general quadrature = 



yields 



{A\k^\) =^[e "0 cos^(/i) + e sin^(/i)^ 



(A7) 
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Thus, one of the quadrature observables (here /i = |) 
gets squeezed, while another quadrature (here fi = 0) 
gets anti-squeezed. Transforming back to the Heisenberg 
picture amounts to (assuming again /? <C 1) adding a 
time dependent phase /i ^ /i + QqI. As a result the 
parametric oscillations lead to squeezing in a quadrature 
which rotates. This explains the oscillations (see Fig. [4]) 
in the fluctuations of the observable we are interested in, 
namely An = A|/c^|, and gives rise to the expected decay 
of the envelope in Eq. (fT4|). 



with k q, and 

i|C(r)) = 7i^|C(r)). (B4) 

When deriving the adjoint equations, we first look at 
simphfied eom and then make the equations subsequently 
more complete. We start with two simple Schrodinger 
equations with projectors, than we add the nonlineari- 
ties, and finally we take into account also the number 
distribution dynamics. 



APPENDIX B: OPTIMALITY SYSTEM FOR 
MCTDHB 

In this appendix we derive for the MCTDHB method 
the control equation, the terminal conditions for the ad- 
joint variables, and the equations of motion (eom) for the 
adjoint variables. The OCT approach is the same as in 
Sec.|IIH where the Frechet derivatives with respect to the 



variables 



2, C, the adjoint variables 



2, C, and 



the control A are equated to zero. 

The Lagrange functional is given analogous to Eq. ([of : 

L(0g,0e,C,0g,0e,C, A) = J(0g,0e,C, A) + (Bl) 

Re{4)k, i4>k - P [Hk + Uq{p}iI ^ Pksqi(t)l(t)q(t)i\ ) 

k s,q,l 

-\-Re{C,iC - {-Q.Jx + ^ XI ^k^\^l^rnWkqlm)C) , 

k,q,l,'m 



where the brackets again imply integrating over time. 
The cost functional J(0g, 0e, C, A) is given in Eq. (23). 
The control equation reads 



VJ = -7A-^i?e( 



^fe|^|0fc> 



(B2) 



+ Re ^ {^q I 0fc > (0fc I ^ I 0q / 
q,k 



-i?e(C|4|C>5^(-l)'' (0.1-^10, 



where (—1)^ = 1 and (—1)^ = —1. In case of formu- 
lation we get a Poisson equation for the gradient, with 
the rhs from above. 

When calculating the Frechet derivatives, partial time 
integration has to be performed, which yields the 
terminal conditions for the adjoint variables. Here 
it is necessary to make the symmetric definition 
d = JJ^dx(i)*g[x\(i)e[x\ - dx(i)*g[x\(i)e[x\)/2, cf. sec- 
tion |lVB2l We then have 



iMT) = -l2{c^i\MT))4 
+7i(C(r)|{^,aIa,}|C(r)> 



(B3) 



Two Schrodinger equations with projectors 

We consider the eom = P/i0^ with Lagrangian (we 
omit arguments for simplicity) 



^ = ^ + ^f 



dtL , with L = Re{(j)i\i^i — Phcpi) . 

(B5) 

We then need to calculate the rhs of = —2-^. With- 
out the projectors we would be left with Schrodinger 
equations for 0^, but here we have in addition the terms 



•i)(Pk 



Two nonlinear Schrodinger equations with projectors 



Here we consider the Eqs. (Bl), but with constant 



coefficients. The Lag rangian is then given by the ffist 
two lines of Eq. (Bl), assuming constant one- and two 



body reduced densities. Without projectors the relevant 
terms (to be multiplied by the corresponding coefficients) 



are 



and 



9lrst = (Plf/^rK^it + (1)1 Sis) • 

The projector terms yield 



(B6) 
(B7) 



Gist = (0d0z>(/>;0.0t + J2(^l\^k)(|)Mt6^r , (B8) 

k 



and 



Irst = 9l\(Pr9s' 



bS^) + Y.^4>k\4>. 



'I (Pr 



'5is + (I>*s5it) . 

(B9) 



MCTDHB equations 



Now we consider the complete MCTDHB equations. 



?(^) ? with Lagrangian given in Eq. (IBII). For the eom for it 
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requires also to consider derivations of the number distri- 
bution part in the Lagrangian, i.e., of the tunnel coupling 
and the two-body matrix elements. This yields the con- 
tributions 

M, = {-iy2Re{C\J,\C)h(f), 

+ Uo/2 {{c\dlaldidm\c) + {c\dldldidm\c)) x 

k,q,l,m 

((Pl^ik + (t)l5iq)(i)l(i)m (BIO) 

Using also the previous results we obtain as eom 



= h(t)i - Si 

+ X] ^^^P^kk {^Pkkkk {glkkk + 9kkkk - Gkkkk - G\ 



+2pkjkj{9kjkj + 9kjkj - Gljkj - G\jkj) (Bll) 

+Pkkjj{9kkjj ~ ^kkjj) + Pkkjji9kkjj ~ ^fc/cjj)} + ■ 

The eom for C is, besides a two-mode Hamiltonian con- 



tribution as in Eq. (20), given by the derivatives of the 



densities in the orbital part of the Lagrangian as 



i\c) = n\c) 



(B12) 



+4i?e(0,|P|0,f (; 



{p^^dld]d,dj 



d]di 



Pijij) 



+2i?e[(0,|P0*0|>(priaIaIa,a, - ^P^.n)]} 



\C). 



APPENDIX C: DETAILS ON THE NUMERICAL 
IMPLEMENTATION 

In this appendix we give details on our numerical im- 
plementation of the time evolution and optimal control of 
the MCTDHB equations. Given the initial value problem 
y = f{yii)i vW = Uoi we employ a Modified Crank- 
Nicolson time-stepping method [52j 



y 



,n+l 



y" 



At 



(CI) 



which is a norm preserving scheme. At each time step, 
we solve this equation using Newton iterations, where 
y"^ converges quadratically and we employ the stopping 



criterion — < 10~^. This time-stepping scheme 

is very stable, even when replacing y'^+^ on the rhs 

(with loss of accuracy off course). 

For Newton's method we have to take functional 
derivatives of our eom, and due to the projector terms we 
obtain equations of type {A+uv^)x = 6, Ax-\-{x'^v)u = 6, 
and Ax + ^{v^x)u = 6, to be solved for the vector x. Here 
^ is a sparse matrix and can thus be inverted efficiently, 
and u and v are vectors. The first equation is linear in x 
and the term uv* , which arises due to the projectors is a 
low rank dense matrix, and its inversion is thus compu- 
tationally very expensive. A way around the direct nu- 
merical inversion is the Sherman-Morrison formula [58] . 
which reads 



x = A-^b- 



v^A-^b 
1 + v'fA-^u 



A 



-1 



(C2) 



and is solvable whenever A is nonsingular and v'^A'^u ^ 
— 1. For computational efficiency, A'"^ is never calculated 
exphcitly, but rather the vectors A~^h and A~^u. 

The second and third equations to be solved are not 
hnear in x due to complex conjugation. We can attack 
this problem by writing every component in terms of real 
and imaginary parts and use the matrix representation 
of complex numbers 



X 7 X nr 



(C3) 



We then have to solve the systems 



Aj, Ai 

Ai Ar 



Xi 



and 



Aj- Ai 

Ai Ar 



b. 
(C4) 



, (C5) 



respectively, by use of the Sherman-Morrison- Woodbury 
(SMW) [58] formula 



{A + UCV)-^ = A-^ - A-^U{C-^ + VA-^U)-WA-^ . 

(C6) 

Here, A^ [/, C and V are matrices of size n x n x 
k X k and k x n, respectively, with k < n. 

The complex extensions to the SMW-formula, which 
are novel, to our knowledge, are also used for the adjoint 
equations (see appendix |A|), which are hnear in the ad- 
joint variable, but again contain full matrices due to the 
projectors in Eqs. (16). 
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